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PREFACE 


This  Memorandum  continues  RAND’s  research  into  the  statistical 
analysis  of  computer  simulation  experiments.  The  overall  purpose  is 
to  find  methods  for  efficiently  extracting  useful  information  from 
time  series  generated  by  these  experiments.  This  particular  study 
describes  a  technique  for  estimating  the  required  sample  size  in  a 
simulation  experiment,  and  provides  flow  charts  and  computer  pro¬ 
grams  for  incorporating  the  proposed  technique  directly  into  a  com¬ 
puter  simulation  program.  Emphasis  is  on  relieving  the  investigator 
of  the  need  to  Interact  with  the  ongoing  simulation  to  determine  when 
the  desired  statistical  precision  has  been  obtained. 

Preceding  work  on  this  subject  is  described  in  G.  S.  Fishman  and 
P.  J.  Kiviat,  Spectral  Analysis  of  Time  Series  Generated  by  Simulation 
Models,  The  RAND  Corporation,  RM-4393-PR,  February  1965;  G.  S.  Fishman 
Problems  in  the  Statistical  Analysis  of  Simulation  Experiments:  The 
Comparison  of  Means  and  the  Length  of  Sample  Records,  The  RAND  Corpora 
tion,  RM-4880-PR,  February  1966;  and  G.  S.  Fishman,  Digital  Computer 
Simulation:  The  Allocation  of  Computer  Time  in  Comparing  Simulation 
Experiments,  The  RAND  Corporation,  RM-5288-1-PR,  October  1967. 
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SUMMARY 

A  method  will  be  described  for  estimating  and  collecting  the 
sample  size  needed  to  evaluate  the  mean  of  a  process  (with  a  spec¬ 
ified  level  of  statistical  accuracy)  in  a  simulation  experiment.  A 
procedure  is  also  described  for  incorporating  the  determination  and 
collection,  of  the  sample  size  into  a  computer  library  routine  that  can 
be  called  by  the  ongoing  simulation  program. 

We  present  the  underlying  probability  model  that  enables  us  to 
denote  the  variance  of  the  sample  mean  as  a  function  of  the  autore¬ 
gressive  representation  of  the  process  under  study.  And  we  describe 
the  estimation  and  testing  of  the  parameters  of  the  autoregressive 
representation  in  a  way  that  can  easily  be  "built  into"  a  computer 
program. 

Several  reliability  criteria  are  discussed  for  use  in  determin¬ 
ing  sample  size.  Since  these  criteria  assume  that  the  variance  of 
the  sample  mean  is  known,  an  adjustment  is  necessary  to  account  for 
the  substitution  of  an  estimate  for  this  variance.  It  is  suggested 
that  Student's  distribution  be  used  as  the  sampling  distribution, 
with  "equivalent  degrees  of  freedom"  determined  by  analogy  with  a 
sequence  of  independent  observations. 

A  bias  adjustment  is  described  that  can  be  applied  to  the  begin¬ 
ning  of  the  collected  data  to  reduce  the  influence  of  initial  condi¬ 
tions  on  events  in  the  experiment.  Four  examples  are  presented  using 
these  techniques,  and  comparisons  are  made  with  known  theoretical 
solutions.  Finally,  we  present  the  minimum  variance  unbiased  esti¬ 
mator  of  the  sample  mean,  which  turns  out  to  be  a  function  of  the 
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autoregressive  coefficients.  Before  these  results  can  be  used  in 
practice,  more  will  have  to  be  known  about  their  sampling  properties. 

In  conclusion,  it  is  noted  that  the  use  of  the  procedures  de¬ 
scribed  here  relieves  the  user  of  the  task  of  continually  interacting 
with  the  simulation  experiment 
within  an  acceptable  range. 


to  determine  whether  his  results  are 
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1 .  INTRODUCTION 

In  sanv  simulation  expe»**eents ,  observations  collected  on  the 
process  of  interest  are  positively  correlated.  This  oeans  that  an 
observation  that  exceeds  (falls  below)  the  population  oean  of  the 
process  tends  to  be  followed  by  another  observation  that  exceeds 
(falls  below)  the  oean.  In  a  sore  liaited  set  of  siculation  experi¬ 
ments,  successive  observations  are  negatively  correlated  so  that  an 
observation  exceeding  (falling  below)  the  oean  tends  to  follow  one 
falling  below  (exceeding)  it. 

Estimating  a  population  aean  fros  scrapie  data  is  a  common  objec¬ 
tive  of  the  statistical  analysis  of  a  simulation  experiment  and,  more¬ 
over,  an  estimate  of  the  variance  of  the  sample  mean  is  helpful  in 
assessing  how  representative  the  saspit  mean  is  of  the  population  mean. 
For  uncorrelated  observations,  the  sample  population  variance  divided 
by  the  sample  size  provides  a  convenient  estimate  of  the  variance  of 
the  saaple  aean.  For  correlated  data,  the  variance  of  the  sample  aean 
is  a  function  of  the  correlation  between  observations,  a  fact  that 
causes  considerable  difficulty  in  estiaating  this  variance. 

Estimators  that  cake  account  of  the  correlation  to  varying  extents 
have  been  suggested  in  the  literature  14. 6j,  but  all  require  a  degree 
of  subjective  judgeent  regarding  their  adequacy.  Ideally,  one  wants 
an  algorithm  that  can  be  "built  into"  a  computer  simulation  program 
and  can  objectively  estimate  the  sample  size  needed  to  obtain  a  speci¬ 
fied  confidence  interval  for  a  population  mean.  Such  a  procedure  would 
relieve  an  investigator  of  the  burden  of  estimating  tbe  variance  of 
the  saaple  aean  frea  a  data  saaple  obtained  from  a  trial  run,  estiasting 


the  sample  size  necessary  for  the  specified  confidence  interval,  and 
then  collecting  that  many  more  observations  in  a  successive  simula¬ 
tion  run.  An  ideal  program  would  accomplish  these  tasks  without  tak¬ 
ing  the  simulation  off  the  computer.  This  Memorandum  describes  an  al¬ 
gorithm  for  doing  this. 

Estimating  the  population,  mean  and  estimating  the  variance  of 
this  resulting  estimate  are  problems  in  statistical  Inference  that 
require  an  underlying  probability  model.  Some  models  are  more  con¬ 
venient  than  others,  and  it  is  natural  to  assume  a  model  that  yields 
desirable  statistical  properties.  For  example,  we  may  assume  that 
observations  are  independent  and  identically  distributed  if  we  be¬ 
lieve  that  the  outcome  of  any  trial  is  not  influenced  by  the  outcomes 
of  other  trials  and  also  if  the  ordering  of  the  trials  does  not  affect 
their  outcomes. 

If  we  suspect  that  the  observations  are  statistically  dependent 
but  that  the  dependence  is  strictly  nonlinear,  then  we  may  assume 
that  the  observations  are  uncorrelated.  This  model  is  less  restric¬ 
tive  than  that  for  the  independent  case,  for  it  implies  that  the  ob¬ 
servations  share  a  coianon  mean  and  a  couaon  variance;  the  covariance 
between  any  two  observations  is  zero,  but  no  specification  is  made 
regarding  the  nonlinear  statistical  relationship  among  observations. 

The  assumptions  of  independent  or  uncorrelated  observations  ap¬ 
ply  in  many  statistical  analyses;  in  simulation  experiments,  however, 
observations  are  often  autocorrelated.  This  means,  statistically, 
that  an  observstion  is  linearly  dependent  on  preceding  observations. 
Since  failure  to  acknowledge  this  autocorrelation  can  seriously  impair 
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the  veracity  of  conclusions  based  on  statistical  results,  a  probabil¬ 
ity  model  is  needed  to  account  for  autocorrelation. 

The  commonly  employed  estimator  of  the  population  mean  for  auto- 
correlated  observations  has  the  same  algebraic  form  as  that  for  inde¬ 
pendent  observations.  In  the  Independent  case,  it  is  statistically 
unbiased  and  has  minimum  variance;  in  the  autocorrelated  case,  these 
properties  hrld  only  for  large  samples.  Nevertheless,  the  estimator 
has  much  intuitive  appeal  and  can  be  computed  easily,  two  properties 
that  account  for  its  common  use. 

The  variance  of  the  sample  mean  is  a  function  of  the  autocorre¬ 
lation  between  observations  and,  therefore,  specification  regarding 
the  covariance  structure  is  necessary  to  make  the  estimation  of  this 
variance  a  tractable  statistical  problem.  The  probability  model  of 
a  covariance  stationary  sequence  provides  a  convenient  framework  with¬ 
in  which  this  problem  and  many  others  can  be  solved.  It  is  this  model 
that  we  describe  in  Sec.  2. 

In  [4],  the  general  properties  of  a  covariance  stationary  se¬ 
quence  enable  us  to  derive  a  useful  estimator  of  the  variance  of  the 
sample  mean;  but,  unfortunately,  that  estimator  is  difficult  to  build 
into  a  simulation  program.  By  sdding  several  mild  restrictions  to  the 
model,  we  may  represent  a  present  observation  as  a  linear  combination 
of  past  observations  plus  a  random  residual  uncorrelated  with  past 
observations.  This  scheme  is  called  an  autoregressive  representation 
of  the  sequence,  and  the  weights  in  the  linear  combination  are  called 
the  autoregressive  coefficients.  For  large  samples,  it  is  shown  that 
a  knowledge  of  the  autoregressive  coefficients  and  the  variance  of 


the  uncorrelated  residuals  enables  vis  to  approximate  the  variance  of 
the  sample  mean  closely. 

While  the  autoregressive  coefficients  and  the  residual  variance 
are  unknown,  we  can  estimate  them  as  described  in  Sec.  3.  We  then  use 
these  to  estimate  the  variance  of  the  sample  mean.  This  approach  is 
desirable  because  the  estimation  and  testing  of  the  coefficients,  the 
computation  of  the  sample  residual  variance,  and  the  estimation  of 
the  variance  of  the  sample  mean  can  all  be  accomplished  directly  in  a 
simulation  program  ind  require  no  user  intervention. 

Section  4  discusses  several  criteria  for  determining  sample  size. 
In  some  experiments  va  may  want  the  confidence  interval  to  have  some 
fixed  absolute  width  around  the  mean.  In  others  we  may  require  the 
width  to  be  a  fixed  percentage  of  the  mean.  To  account  for  the  use 
of  an  estimate  for  the  variance  of  the  sample  mean.  Sec.  5  introduces 
the  t  distribution  with  appropriate  adjustments  for  its  use  with  auto- 
correlated  data.  Initial  conditions  often  influence  the  behavior  of 
the  process  under  study;  and  it  is  desirable,  whenever  possible,  to 
reduce  the  extent  of  this  influence.  Section  6  covers  this  topic. 

In  Sec.  7,  several  examples  with  known  solutions  are  presented 
to  illustrate  how  well  the  techniques  work  statistically.  The  exam¬ 
ples  include  zero-,  first-,  and  second-order  autoregressive  schemes 
and  a  single-server  queueing  problem  with  independent  and  exponential¬ 
ly  distributed  interarrival  times  and  service  times.  These  examples 
are  simulated  and  the  estimated  sample  sizes  compared  with  known  theo¬ 
retical  results  presented  in  [5]. 

Earlier  we  remarked  that  the  conventional  estimator  of  the  mean 
of  an  autocorrelated  sequence  la  unbiased  and  minimum  variance  only 


i 


c 

for  large  samples.  It  is  therefore  instructive  to  study  what  can  be 
done  to  derive  an  improved  estimator  of  the  mean  for  a  moderate  sample 
size  and  how  feasible  and  worthwhile  it  is  to  use  the  improved  esti¬ 
mator.  Section  8  presents  the  minimum  variance  estimator  and  compares 
it  with  the  conventional  mean  estimator. 

We  conclude  that  the  algorithms  suggested  here  can  contribute 
significantly  to  solving  the  problem  of  determining  sample  size  in  a 
simulation  experiment.  Because  they  can  be  used  while  minimally  in¬ 
volving  the  simulation  experiment  itself,  they  are  worthy  of  consider¬ 
ation,  especially  since  they  can  be  easily  modified  to  meet  individ¬ 
ual  user  preferences.  Also,  the  suggestions  regarding  reliability 
criteria,  unbiasedness,  minimum  variance  and  variance  reduction  tech¬ 
niques  provide  users  with  information  that  enables  them  to  draw  use¬ 
ful  inference  about  the  process  being  studied. 
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2.  THE  MODEL 

In  many  simulation  experiments,  the  process  of  interest  appears 
as  a  sequence  of  events  in  which  the  index  that  orders  the  events  may 
play  a  role  in  defining  the  relationship  among  events.  The  index  may 
be  time;  for  example,  the  number  of  units  waiting  for  service  at  time 
t,  where  t  assumes  T  different  values.  The  index  may  simply  denote 
order;  for  example,  the  waiting  tine  for  the  t**1  unit  to  receive  ser¬ 
vice,  where  t  assumes  the  values  1,  . ...  N. 

When  the  set  of  ordered  events  is  subject  to  random  variation, 
it  is  called  a  stochastic  sequence.  Let  X£  be  the  value  assumed  by 
the  tth  event  in  which  the  index  t  runs  over  the  Integers.  Then  we 
denote  the  stochastic  sequence  by  {Xt;  t  -  0,  +  1,  . ..,  ±  •}  or,  mere 
concisely,  by  X.  We  could  index  the  events  on  a  finite  set  of  non¬ 
negative  Integers,  but  the  above  definition  of  X  offers  several  expo- 
sltional  conveniences  without  impairing  its  applicability  in  the  pres 
ent  context. 

Let  the  sequence  X  have  mean 

(2.1)  -  E(Xt) 
and  autocovariance  function 

(2.2)  Ra  t  -  E[(Xs  -  ^)(Xt  -  . 

If  is  finite  and  independent  of  the  ordering  index  t,  and  R  is 
c  s,t 

finite  and  a  function  only  of  the  number  of  intervening  events  s-t, 
then  we  may  write  the  mean  and  autocovariance  function  as 
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(2.3) 


*  =  , 

Rs-t  -  Rt-s  =  Et(Xs  '  *><Xt  -  M  * 


respectively.  A  sequence  satisfying  (2.3)  is  called  covariance,  mean- 

* 

square,  weakly  or  wide-sense  stationary. 

Suppose  that 


s  *  t 

s  4  t 

Then  X  is  a  sequence  of  uncorrelated  events,  and  conventional  methods 
of  statistical  inference  apply  when  estimating  p.  In  general,  (2.4) 
does  not  hold  and  more  sophisticated  inferential  methods  are  needed. 

For  a  sample  of  N  observations,  we  compute  the  conventional  sample 
mean  as 


(»2 

t° 


N 

(2.5a)  » 

t*l 


See  [l)  for  a  more  complete  description  of  covariance  stationary 
sequences . 


with  variance 


N-l 

Var^)  =  (1/N2)  ^  Rst 

s ,  c=l 


(2.5b) 


N-l 

=  (1/N)  ^(1  -  |s|/N)Rs  . 

s=l-N 


He  also  require  that 


(2-6) 


lig  R  =  0  . 
s-3®  s 


This  restriction  is  reasonable,  for  we  would  expect  the  covariance 
between  events  in  the  series  to  vanish  as  the  number  of  intervening 
events  increases.  Then  one  may  show  that 


(2.7a) 


P2 


N-l 

IX 


s»l-N 


m  <  •  » 


N-l 

(2.7b)  ££g  ^  (Uj/N)Rs  =  0  , 

s=*l-N 


so  chat  for  large  N 
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Condition  (2.6)  requires  the  absence  of  any  regularly  periodic  com¬ 
ponents  in  X.  If  any  were  present ,  then  the  autocovariance  function  R 
would  contain  undamped  cosine  terms  that  would  violate  (2.6)  and  prevent 
the  convergence  of  (2.7b) .  The  following  example  demonstrates  the  truth 
of  this  assertion.  Let  Xt  be  defined  by 


(2.9) 


X^  =»  a  sin  bt  +  y  , 


where  a  is  a  random  variable  with  zero  mean  and  unit  variance ,  b  is  a 
constant,  and  y  is  a  covariance  stationary  sequence  with  mean  v  and 
autocovariance  function  P.  Moreover, 


(2-10) 


lim  P 
s-w  s 


0. 


Then  X  has  the  autocovariance  function 


(2.11) 


R  =  cos  bs  +  P_ 

S  s 


and,  for  a  sample  of  N  observations. 


N-l 

(2.12)  VaKX^  -  (1/S )  yM  (1  -  |  si /N) (cos  bs  +  Pg)  - 

s=l-N 

Now 


N-l 

cos  bs  »  sin  ib(2N  -  l)/2]/sin  (b/2)  , 

s=l-N 
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N-l 

^  *  Jsjcos  bs  =*  N  sin[b(2N-l)/23/sin(b/2)  -  [l  -  cos(Nb)3/[2sin^(b/2) ]> 
s-l-N 


so  that 


N-l 

(1  -  jsJ/N)  cos  bs  ■  [1  -  coa  (bN)]/[2N  sin^  (b/2)]  . 

s-l-H 


Then 


(2.13)  Jig  VarCiy  -  ^  )  £*,  , 

so  that  (2.8)  does  not  hold  in  this  case.  This  result  suggests  that 
any  regularly  periodic  components  in  the  process  X  be  removed  prior 
to  estimating  the  mean  u- 

The  reader  may  wonder  why  (2.8)  is  of  such  great  significance, 
since  m  is  an  infinite  sum  of  autocovariances.  The  answer  is  that  m 
can  be  computed  from  alternative  formulae  wherein  the  individual  ' s 

need  not  be  known.  As  a  result  our  emphasis  is  on  the  quantity  m  and 
its  estimation.  The  autoregressive  representation,  which  is  to  be  intro¬ 
duced  shortly,  provides  alternative  formulae  for  computing  a. 

Condition  (2.6)  implies  that  the  correlation  between  two  events 

X  and  X  goes  to  zero  as  the  interval  Is-tj  becomes  large.  If  wc 
s  t 

impose  the  added,  but  mild,  restriction  that  there  exists  a  finite 


intervax  r  such  that  two  events,  and  X^,  ar^  statistically  indepen- 
dent  if  js  -  tj  >  r,  and  that  E( | X^_  j  )  <  ®,  then  it  can  be  shown  that 
the  Uniting  distribution  of  the  quantity  tf5  (X^  -  u)  is  noriaal  with 
mean  zero  and  variance  m  [7,  pp.  215-219].  Hereafter,  we  assume  that 
S  is  sufficiently  large  for  us  to  use  this  Halting  result. 

One  of  the  desirable  features  of  a  covariance  stationary  sequence 
is  its  connection  with  a  sequence  of  uncorrelated,  identically  distri¬ 
buted  random  variables.  Using  the  Wold  decomposition  theorem  {16],  one 
say  write  for  X^  satisfying  (2.3)  and  (2.6)' 


(2.14) 


X.  »  a  +  >  a  X 
£  /.  s  t-s  * 


where  ia^ ;  s=0,+l,+2, —  is  a  sequence  of  real  constants  with 


(2.15) 


2  s <  ■ : 


and  (Y^;  t  *0,  t  1,  t  2,  ...,  ±  *}  is  a  sequence  of  uncorrelated,  iden 

2 

tically  distributed  random  variables  with  mean  zero  and  variance  c  .  T 
is  an  appealing  form,  for  we  can  now  write  the  autocovariance  of  X  as 


(2.16) 


R  =  o  7  ,  a  a 
s  f  4  t  s+t  * 


For  a  concise  description  of  the  Wold  decomposition  theorem, 
see  [2,  pp.  286-2881. 
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b  *  0  s  <  0  » 
s 

(2.23) 

Lb  X'  *  Y  . 
s  t-s  t 

3*0 

Moreover,  the  movitig  average  (2.14)  is  expressible  as 

00 

<2-24>  K  “XkVs  • 

s*0 

Expression  (2.24)  has  an  intuitive  appeal,  for  it  implies  that  is 
a  linear  combination  of  uncorrelat.  d.  identically  distributed  present 
and  poet  events. 

The  number  of  coefficients  in  the  b  sequence  remains  to  be  con¬ 
sidered.  Suppose  that 

bs=0  g  >  p  >  0  , 

so  that 

P 

<2-25>  a  '  • 

3*0 

One  would  normally  expect  thct  after  some  lag  p,  the  contribution  to 
the  behavior  of  Xt  made  by  variables  X  _j,  X  _2,  •••  would  be 
negligible,  so  that  (2.25)  would  adequately  describe  the  relationship 
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betweer.  X'  and  Y.  This  means  that  we  can  find  a  linear  combination 

of  present  and  past  values  of  X'  that  form  a  sequence  Y  of  uncorrelated* 

2 

identically  distributed  events  with  mean  zero  and  variance  a  . 

Using  (2.25)  leads  to 


2,.  2 

m  »  a  /b  , 

(2.26) 

P 

b  -£>s  . 

s=Q 


To  estimate  a  and  the  coefficients  in  the  b  sequence,  we  apply  the 

linear  least-squares  method  to  a  sample  of  observations  on  X  using  the 

autoregressive  representation  (2.25).  We  can  subsequently  estimate 

2 

m  by  substituting  estimates  of  o  and  the  b  *s  into  (2.26).  The  sample 

A  ^ 

variance  of  the  sample  mean  is  then  m/N. 

It  is  instructive  to  compare  m  with  N  Var (X^)  to  measure  the  ade¬ 
quacy  of  the  approximation.  Suppose  X  has  the  autoregressive  represen¬ 
tation 


(2.27)  Xt  -  gX^  -  Yt  US  <  1  , 


so  that 
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1  Rg  =  cr2g^/(l-g2) 

(2.28)  /  H  VarCy  =  l/(l-g)2  -  2g(l-gN)/rN(l+g)(l-g)3] 

m  =  l/(l-g)2  . 

Notice  that  the  restriction  |gj  <  1  is  equivalent  to  requiring 

(2.29)  B(z)  =*  1  -  gz  a  0 


to  have  its  root  outside  of  the  unit  circle. 

Figure  1  shows  the  ratio 

(2.30)  q  -  [Vat(XN)/VN]% 

for  several  values  of  g  and  varying  sample  sizes  N.  The  square  root 
comparison  is  appropriate,  for  it  is  the  standard  deviation  of  the 
sample  mean  that  determines  the  width  of  a  confidence  interval  for 
the  population  mean.  Notice  that  for  Jgj  <  0.50,  the  error  of  approxi¬ 
mation  is  less  than  5  percent  for  N  >  25.  For  Jg|  =  0.95,  the  error 
is  about  10  percent  for  N  >•  100. 

The  error  patterns  differ  noticeably  for  positive  and  negative 
values  of  g.  For  positive  values,  a  always  overestimates  N  Var (X^); 
for  negative  values,  it  always  underestiaates  it.  Also,  N  Var  (jy 
oscillates  when  g  <  0,  the  pattern  being  most  apparent  for  small  N. 

From  inspection  of  Fig.  1,  it  is  clear  that  choosing  an  even  value 
of  N  improves  the  approximation. 


Fig.  1 — Comparison  of  sample  mean  variance  and  approximation  for  a  Markov  process 
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Suppose  we  write 


Jb  X' 

S  t-3 


b  K 

S  t- 


+  Y 

3  t 


3=0 


s=r+l 


(2.31) 


b0  "  1;bs  “  °»  a  >  p  ’ 


so  that 


(2.32) 


E(X '  ,Z) 

t-r-1  t 


Eh  E(X*  X'  .) 
s  v  t— s  t-r-r 

s=0 


VhR 

/  j  s  s-r-i 
s-0 


Here  Zfc  is  a  conditional  random  variable  derived  by  removing  the 
linear  effects  of  X£  x^-r  *ro®  Xt'  r  3  P»  then  all  the  linear 

effects  of  past  events  have  been  removed  from  X£  so  that 


(2.33) 


E<Xt-p-lZt>  "  0  • 


b  *  0  • 
P+-1 


This  result  is  important,  for  it  gives  ua  a  way  of  determining 
p,  the  order  of  the  autoregressive  representation.  Suppose  we  esti¬ 
mate  the  coefficients  in  R  autoregressive  schemes  of  orders  1,  . R. 
He  then  test  the  significance  of  the  last  estimated  coefficient  in 
each  scheme.  By  choosing  R  sufficiently  large,  we  can  find  a  value 
p  <  R  beyond  which  all  remaining  last  estimates  are  insignificant. 


-19- 


It  is  important  to  check  a  number  of  schemes  for  significance 
rather  than  proceeding  stepwise,  successively  testing  the  last  esti¬ 
mate,  and  stopping  when  the  coefficient  is  insignificant.  It  nay 
occur  that  a  coefficient  br  for  R  <  p  vanishes.  Proceeding  in  a  step 
wise  fashion  would  cause  us  to  stop  testing  prematurely. 
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3 .  STATISTICAL  INFERENCE 

Consider  a  set  of  N  observations,  X^,  ....  X^,  that  we  wish  to  fit 
to  an  autoregressive  representation: 

r 

(3.1)  •  *t  • 

s-0 

We  first  adjust  for  the  sample  mean, 

0.2)  x;  -  xt  -  \  * 
and  then  compute  the  sample  autocovariances, 

N-t 

(3.3)  c^  r  -  (1/N)  y\;x;^  t  -  o,  l,  r. 

t-i 

The  quantity  R  is  the  highest  order  of  the  autoregressive  schemes  that 
we  plan  to  test.  Since  we  plan  to  compute  estimates  for  schemes  of 
successively  higher  order,  we  use  a  tine -conserving ,  recursive  estima¬ 
tion  procedure  suggested  by  Durbin  [3],  which  Whittle  also  describes 
[15,  p.  37}. 

A  ^ 

Let  b  .  be  the  s  lagged  estimated  coefficient  in  the  scheme 
r+i  |S 

of  order  r+1,  and  let 

A 

br+l,0  " 


(3.4) 


1 


r  ■  0 ,  1, 


R-l. 


Then  from  the  r  order  scheme,  where  r  =  0,  1,  . R-l ,  we  compute 


,  .y:  c 

r  r,s  N,j 


'r  -E"br,sC»,r- 


fa  ,  -  ,,  =  -w  Iv 

r+1 ,  r+1  r  r 


fa  .  =  b  +  b  ,,  .,b  s  =  l,  — ,r 

r+i,s  r,s  r+1, r+1  r,r-s+l 


The  sample  residual  variance  that  is  an  estimate  of  o  is 


(3-3) 


N  r  r+1 

o2  =  (B  -  r-1)-1^  ^ fa  (X 

r+1  '  /  /  /  -<  r+1  ,s  t-s 

t=r+l  *-s=0 


(3-6) 


vr  r  =  0’  — • R- 1 


Using  a  result  based  on  Whittle  fl5,  pp.  72-731,  one  may  show  that 


for  large  N 


var(b  >  ~  yu  /N 
r,r  r 


w  =»  1  -  fa 
r  r,r 


If  the  confidence  interval 


(3.7) 


bc,r±F(»r/f0: 


covers  zero,  then  we  accept  the  hypothesis  that,  fa  =  0.  The  quantity 
P  is  the  point  on  the  normal  curve  corresponding  to  a  significance  level 


cr ,  where 
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(3.8) 


1  -  a/2  . 


This  method  of  determining  p  is  based  on  the  testing  procedure 
described  by  Jenkins  and  Watts  in  [8,pp.  189-2003-  There  the  test 
statistic  has  the  Student  t  distribution  but,  in  the  present  context, 
we  assume  N  to  be  sufficiently  large  so  that  the  normal  approximation 
is  acceptable.  Quecoullle  [13]  has  described  an  alternative  and  sore 
precise  large  sample  :,goodneaa  of  fit"  teat  for  autoregressive  schemes 
Unfortunately,  his  test  requires  several  more  complicated  computations 
than  thoae  described  above,  and  it  does  not  appear  to  ba  easily  in¬ 
corporated  into  a  simulation  program. 

Suppose  thst  b  is  tested  for  significance  for  r  «  1,  ...,  R, 
r  ,r 

and  that  the  coefficient  b  ,  p  <  R  is  significant,  but  the  coeffi- 

P»P 

dents  b  are  not  significant  for  r  •  p  +  1,  ...,  R.  Then  we  choose 
the  order  of  the  scheme  to  be  p  and  estloate  a  by 


(3.9) 


Three  questions  remain  to  be  answered  here.  One  is  the  choice 


of  R;  the  second,  the  choice  of  the  initial  sample  size;  and  the  third 
the  choice  of  a  and  thereby  P.  The  larger  is  R,  the  better  is  the 
chance  of  including  the  correct  p  within  the  schemes  tested.  But  low- 
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order  autoregressive  schemes  often  suffice  to  account  for  the  auto¬ 
correlation  structure  ir  X  and,  consequently,  choosing  R  to  be  10  should 
be  more  than  adequate.  In  our  work  we  have  chosen  R  to  be  4  and  5. 

The  steps  that  we  have  so  far  described  are  based  on  an  initial 
sample  of  N  observations.  For  the  normal  approximation  to  be  acceptable* 
we  require  that  N  -  R  >  30  so  that  the  initial  sample  size  il  exceeds 
R  +  30;  how  much  greater  it  should  be  depends  on  the  cost  of  collecting 
observations,  a  point  to  which  we  return  in  Sec.  6. 

Earlier  we  spoke  of  testing  the  last  coefficient  in  each  of  the  R 
autoregressive  schemes  estimated.  We  then  have  a  multiple  testing  prob¬ 
lem  and,  if  we  choose  a  significance  level  a  for  each  test,  the  signi¬ 
ficance  level  for  the  multiple  test  will  be  greater  than  a.  The  greater 
R  is,  the  greater  the  divergence  is  between  a  and  the  significance  level 
for  the  multiple  test.  It  is  therefore  recommended  that  the  choice  of 
or  be  less  than  one  would  customarily  use  in  testing  a  single  hypothesis. 
This  divergence  is  also  a  good  reason  for  keeping  R  snail. 
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4.  RELIABILITY  CRITERIA 


The  purpose  a£  Che  presenc  research  is  to  deteraine  Che  sample 
size  needed  Co  estimate  Che  populaCion  mean  with  a  specified  accuracy 
or  reliabilicy.  Since  we  are  creacing  Che  sample  mean  X^,  as  a  normal 
variaCe  wich  approximaCe  variance  a/ll,  we  specify  reliabilicy  by  mean  £ 
of  the  confidence  stacemenc 

(4*1)  Pr^*N  "  *1  <  i  -  3  » 

where  the  confidence  level  3  is  a  saall  probability  such  as  0.05  or 
0.10,  and  Q  is  the  normal  point  corresponding  to 


(2tt) 


-1/2 


We  say  also  write  (4.1)  as 

(4.2)  Pr(Xj.  -  <  P  <  \  +  Q^/N)  ~  1  -  3  , 

and  we  note  that  the  larger  is  N,  the  shorter  is  the  confidence  interval 
around  u. 

Suppose  we  wish  to  collect  a  sample  size  such  that  the  variance  of 
the  resulting  sample  mean  is  less  than  or  equal  to  V  with  probability 
1-9.  That  is, 

(4.3)  Pl(  <  u  <  V  +  ~  -  3  • 

To  determine  K*,  we  note  the  equivalence  of  (4.2)  and  (4.3)  so  that 
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(4.4) 

V  =*  s/ff* 

» 

(4.5) 

N*  =»  a/V  . 

If  c  were  known  a  priori,  then  the  determination  of  N*  in  (4.5) 
weald  follow  directly.  Since  we  only  have  an  estimate  fi  ->£  m,  it  is 
natural  to  replace  m  by  n  in  (4.5)  and  so  estimate  N*.  As  a  first 
approximation,  this  is  a  reasonable  approach.  Figure  2  shows  a  flow 
chart  oi  an  iterative  procedure  for  determining  and  collecting  K*  ob¬ 
servations. 

Notice  that  we  collect  y(N*  -  S) ,  not  S*  -  N  additional  observa¬ 
tions.  This  scaling  increases  the  number  of  iterations,  but,  more 
important,  it  decreases  the  total  number  of  unnecessary  observations 
collected.  Since  the  computer  time  saving  will  generally  be  such 
greater  with  regard  to  the  avoidance  of  collecting  unnecessary  obser¬ 
vations  than  the  time  expended  on  additional  Iterations,  scaling  is 
desirable.  In  the  examples  to  be  presented  in  Sec.  8,  y  was  set  equal 
to  1/3  and  1/2  for  comparison. 

The  specification  of  V  is  a  statistically  oriented  constraint  in 
terns  of  population  parameters.  Often  we  prefer  a  confidence  statement 

(4.6)  Pr<*K*  -  c<n<X^  +  c)~l-3, 

where  c  is  a  specified  constant.  Here  we  wish  to  determine  a  sample 
size  N*  such  that  the  probability  is  about  1-3  that  the  difference 
between  the  sample  and  population  means  does  net  exceed  ±  c.  This 
is  an  absolute  reliability  criterion  to  be  net.  To  determine  N*,  we 
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note  that  (4-,2)  and  (4.6)  are  identical  for 


so  that 

(4.8)  N*  «  m(Q/c)2  . 

And  we  note  that  the  relationship  between  V  and  c  is 

V  -  (c/Q)2  . 

As  an  alternative,  we  cay  wish  to  determine  N*  such  that 

(4.9)  -  cji  <  p  <  +  cn)  ~  1  -  g  , 

so  that  the  probability  is  1  -  B  that  the  difference  between  the  sample 
and  population  means  does  not  exceed  +  cu.  Using  (4.2)  again,  we  have 

(4.10)  N*  -  m(Q/cv02  ♦ 

This  is  a  relative  reliability  criterion.  Notice  in  (4.8)  and  (4.10) 
that  halving  c  causes  a  fourfold  increase  in  N*.  Here  m  and  ja  are  un- 

A 

known,  so  we  replace  them  In  (4.10)  by  m  and  . 

Both  the  absolute  and  relative  criteria  are  being  determined  in 
(4.8)  and  (4.10)  using  Q  from  the  normal  distribution.  In  the  absolute 
case  we  know  from  theoretical  considerations  that  failure  to  account 

A 

for  the  substitution  of  m  for  tn  makes  Q  smaller  than  it  should  be  and, 
therefore,  causes  an  underestimate  of  N*.  For  the  relative  case,  the 
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problem  !.s  further  compounded  by  the  substitution  of  for  u. 

By  drawing  several  analogies  to  the  case  of  independent  observa¬ 
tions  with  unknown  mean  and  variance,  we  can  introduce  a  correction 
factor  for  Q  to  account  for  the  unknown  m.  This  is  done  in  the  next 
section. 


■HWSBgg  ggg  gwgggg 
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5.  IMPROVED  CONFIDENCE  INTERVALS 


The  mean  and  variance  of  the  process  under  study  are  \i  and  R^»  re¬ 
spectively.  Let  us  new  consider  a  hypothetical  process  made  up  of  in¬ 
dependent  events  each  with  mean  and  variance  R^.  For  L  observations 
the  variance  of  the  sample  mean  of  the  hypothetical  process  is  R^/L. 

In  the  process  under  study  it  is  m/N  for  large  N.  Equating  these  sample 
mean  variances ,  we  have 

(5.1)  Rq/L  =  m/N  , 

so  that,  with  regard  to  the  sample  mean  variance, 

(5.2)  K  *  N/L  -  m/Rjj 

is  the  number  of  observations  to  be  collected  on  the  process  under  study 
that  is  equivalent  to  collecting  one  independent  observation  on  the  hy¬ 
pothetical  process. 

Suppose  we  have  an  estimate  of  R^  for  the  hypothetical  process. 

Then  we  use  the  t  distribution  with  L  -  1  degrees  of  freedom,  together 
with  X^  and  Cq/L,  to  compute  a  confidence  interval  for  u.  The  use  of 
L  -  1  Instead  of  L  is  due  to  the  substitution  of  for  t*  in  C^.  In 
the  present  study  we  may  derive  a  more  representative  confidence  interval 
for  u  by  using  the  t  distribution  with  X  ,  and  C_/L  with  L  -  1  equiva- 
lent  degrees  ol  freedom,  where 


(5.3) 


L  -  N/K  »  NRq/tii  , 
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the  loss  of  one  degree  of  freedom  being  for  t-he  sample  mean  substitution. 

A 

To  estimate  L,  we  replace,  a  in  (5.3)  by  m. 

To  incorporate  a  table  of  the  appropriate  critical  values  for  the 
t  distribution  into  a  computer  program  seems  undesirable  since  this 
would  require  a  value  for  each  number  of  degrees  of  freedom  for  each 
9.  Instead,  we  may  use  formulae  for  the  asymptotic  expansion  of  the 
critical  value  of  t  around  the  critical  value  of  the  normal  distribu¬ 
tion  tor  a  given  0.  The  interesting  characteristic  of  the  asymptotic 
expansion  is  that  it  is  a  power  series  in  inverse  powers  of  the  number 
of  degrees  of  freedom  so  that  we  may  compute  Q  for  a  given  3  simply  by 
inserting  the  number  of  degrees  of  freedom  in  the  formulae.  These  for¬ 
mulae  may  be  found  in  [14,  p.  948]. 

It  is  to  be  noted  that  the  use  of  the  t  distribution  corrects  for 
unknown  m.  At  the  present  stage  of  research  no  correction  can  be  offered 
for  the  unknown  p  in  the  relative  reliability  case.  Kor  is  their  any 
adjustment  for  the  use  of  an  estimate  for  L.  To  check  on  the  extent  of 
degradation  due  to  these  omissions,  our  examples  are  based  on  the  rela¬ 
tive  reliability  criterion. 
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6.  BIAS  ADJUSTMENT 

AC  che  beginning  of  Che  simulation  experiment,  a  number  of  vical 
variables  are  assigned  predece mined  values  to  "prime"  the  system.  This 
procedure  establishes  a  set  of  initial  conditions  and,  because  of  the  in¬ 
herent  dependence  among  events  ,  the  first  observation  on  the  process 
of  interest  is  a  function  of  these  initicl  conditions.  The  second  ob¬ 
servation  X2  is  also  a  function  of  these  values  but  to  a  lesser  extent 
than  is.  Successive  observations  are  less  dependent  on  the  initial 
conditions  so  that  eventually  events  in  the  simulation  experiment  are 
independent  of  them. 

Because  of  their  dependence  on  the  initial  conditions,  observations 
near  the  beginning  of  the  experiment  are  not  representative  of  the  pro¬ 
cess  of  interest  and  their  inclusion  in  X^  makes  this  quantity  a  biased 
estimator  of  the  true  mean  p.  As  N  becomes  large,  the  bias  goes  to  zero 
since  early  observations  become  less  influential  on  the  average.  But 
for  moderate  N,  the  bias  may  be  significant  and  should  be  reduced  if 
possible. 

We  noted  that  K  in  (5.2)  essentially  measures  the  number  of  auto- 
correlated  observations  per  independent  observation.  Intuitively,  we 
therefore  expect  the  correlation  between  observations  X  units  apart  to 
be  low.  For  example,  in  the  first-order  process  described  earlier  we  have 

C6*l)  X  -  (1  +  g)/(l  -  g)  , 

80  that  the  correlation  between  the  first  and  X  +  1st  observations  is 
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(6.2)  g(l+g)/(l  g)  <  l/e2  =  .13534  . 

As  a  first  step  toward  reducing  the  influence  of  initial  condir.ions> 
we  remove  the  first  K  observations  from  the  sample.  Then  we  have  the 
sample  mean 


(6.3) 


N 


t-K+1 


with 

(6.4)  VN  K  «  m/(N-K)  =*  a2/[(N-K)b2 j  =  J  . 

The  confidence  interval  for  u  is  then  computed  using  the  t  distribution 
with  L  -  2  "equivalent"  degrees  of  freedom  and  Xjj  R  and  CQ/(L-1)  ,  L 
being  estimated  as  before. 

Figure  3  is  a  flow  chart  that  illustrates  one  way  of  including 
the  bias  adjustment.  When  the  sample  size  is  judged  sufficient,  the 
sr.-iple  mean  is  recomputed  using  the  newly  •.stimaced  bias  adjustment  K. 

A  A 

The  estimates  m  and  L  are  not  recomputed,  since  our  experience  has  shown 
that  a  recomputation  of  these  quantities  makes  little  difference. 


5«  Collect 

S  r A  /?  r  N  ■  M  L**4  observitlons 


Fig.  3  -  Computation  of  unweighted  sample  mean  with  bias  adjustment  (M,  P,  R,  V  and  Y  specified) 
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7 .  EXAMPLES 

This  section  presents  four  exaaples  to  illustrate  how  the  pro¬ 
posed  technique  works.  The  first  three  examples  are  zero-,  first- 
and  second-order  autoregressive  schemes  for  normal  stochastic  sequences. 
One-hundred  replications  were  collected  on  each  example  to  enable  mean¬ 
ingful  comparisons  for  different  significance  levels  a  in  determining 
the  order  p  of  tne  autoregressive  schemes  and  for  different  values  of 
Y,  the  scale  factor  used  in  collecting  additional  observations.  The 
fourth  example  is  a  first-come,  first-served  single-server  queueing  prob¬ 
lem  with  independent  ansi  exponentially  distributed  interarrival  and 
service  times.  The  purpose  of  this  example  is  to  study  a  problem  more 
closely  akin  lx  chose  usually  analyzed  in  discrete  event  simulation 
experiments  and  also  one  in  which  the  underlying  distributions  are  not 
normal.  Analytical  solutions  are  available  for  all  four  exaaples  and 
serve  as  a  check  on  the  technique.  The  computer  program  was  written 
in  the  SIHSCR1PT  II  programming  language  [10] . 

In  Che  first  example  we  considered  a  stochastic  sequence  where 

(7.1)  Xt  =  0.5+Yt, 

Y  being  normal  with 


(7.2) 


t  ^  s  . 


(7.3) 


S<Xt)  -  *  »  .5  , 


and  for  a  staple  of  N  observations. 


(7.4) 


VarfXjj)  *  1/N 


In  the  second  exaaple  we  considered  a  first-order  scheae 


(7.5) 


Xt  *  °-5Xt-l  +  5  +  yt  . 


so  that 


(7.6) 


E(xt)  »  *  •  1  , 


VarCX^)  ~  4/N  . 

And  for  the  third  exaaple  we  studied  a  second-order  scheae. 


xt  "  0  5Xt-i  '  -25V?  +  0-5  +  r  , 


(7.7) 


H(\)  -  -  2/3 


Var ~  I6/(9N>  . 


”*  ,U*’’UU"  V*r(V  “”*«■>  uslne  (2.26) .  n,  co.fttci.dU  ta 

th.  second -order  sc  hero  chosen  to  i  Host  re  re  ho.  .  hlgh.r-ord.r 

autoregressive  representation  does  not  necessarily  iaply 
correlation  in  the  sequence  and,  hence,  a  larger  Var(X^ . 


aore  auto- 
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The  objective  in  all  three  exaaples  was  to  obtain  sample  sizes 

* 

N  such  that 


PrClx^  -  uj  <  cm)  ~  1  -  3  . 

(7.8)  c  =  0.20  , 

P  =0.10  . 

Using  (4.5)  ,  the  required  sanple  size  for  each  exanple  was 

(7.9)  N*  =  (16. 4)2  ~  269  , 

so  that  the  respective  sample  oean  variances  were 

Var( =  0.003717  , 

(7.10)  Va^X^  ~  0.014870  , 

Var(X^)  -  0.006609  . 

For  each  example ,  two  significance  levels  for  determining  the  auto¬ 
regressive  order  were  studied.  They  were  a  =  0.023  and  a  =  0.05.  Also 
two  scaling  factors  were  examined  for  determining  the  number  of  addi¬ 
tional  observations  to  be  collected.  They  are  v  =  0.3  and  y  =  0.3333. 
Therefore  each  example  contained  four  cases. 

The  results  for  the  total  of  twelve  cases,  each  containing  IOC- 
independent  replications,  are  presented  in  Table  I,  where 
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p  =  order  c-f  the  autoregressive  scheme  (Sec.  2) 
a  =  significance  point  in  test  to  determine  p  (Sec.  3) 
v  =  weighting  factor  for  newly  computed  required  sample  size  (aec.  5} 
VM  =  approximation  to  variance  of  sample  mean  (Sec.  2)  . 

n 

In  all  experiments,  V  was  close  to  Var(X^.). 

Column  4  lists  the  average  required  s.'.mple  size  computed  on  the 
last  iteration  of  each  experiment.  Column  3  lists  the  average  sample 
si2*  collected  on  the  last  iteration,  which  is  naturally  greater  than 
the  corresponding  quantity  in  col.  4  since  it  is  precisely  this  condi¬ 
tion  that  terminates  the  experiment.  The  quantities  n  parentheses  are 
the  saaple  standard  deviations.  The  highest  order  autoregressive  scheme 
R  considered  was  4. 

Notice  that  increasing  a  from  0.023  to  0.05  causes  slightly  less 
than  a  doubling  in  col.  8  for  p  =  0.  A  less  narked  increase  occurs 
for  p  »  1,  2.  Also  noteworthy  is  the  general  increase  in  col.  8  for 
a  given  a  as  p  increases.  These  increases  would  be,  in  fact,  larger 
if  R  were  greater  since  sore  tests  would  be  performed.  Since  we  ex¬ 
pect  p  >  0,  it  is  advisable  to  sake  a  snail  and  also  tc  restrict  R. 

The  choices  of  a  -  0.025  and  R  5  4  appear  to  be  acceptable  operating 
conditions  in  the  cases  described. 

When  y  is  reduced  from  0.5  to  0.333,  the  excess  saaple  size  (col.  6) 
becomes  smaller.  Compared  with  the  theoretical  l«*  of  269,  the  use  of 
a  =  0.5  causes  an  average  of  1.1746  observations  to  be  collected  for 
every  required  observation,  whereas  the  use  of  a  *  0.3333  requires  an 
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Table  1 

TEST  RES l ITS  FOR  100  REPLICATION  5 
(Tieore.ical  N  *269) 
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average  of  1.0195  observations.  The  sample  standard  deviations  in 
col.  5  are  also  smaller  for  y  *  0.3233,  but  the  average  number  of 
iterations  (col.  7)  is  notably  increased.  Moreover,  27  percent  more 
CPU  time  is  required  for  y  *  0.3333.  These  facts  suggest  that  smaller 
y's  improve  statistical  precision  but  require  more  CPU  time.  We  now 
study  the  cost  of  this  improved  precision. 

The  actual  generation  of  data  used  in  these  experiments  consumed 
relatively  little  CPU  time  so  that  we  may  reasonably  attribute  the 
total  CPU  time  to  the  proposed  statistical  technique.  For  y  ■  0.5, 
the  program  processed  274  observations  per  second;  for  y  ■  0.3333,  it 
processed  187  observations  per  second.  In  all  experiments  sumnarized 
in  Table  1,  cha  data  were  retained  in  the  magnetic  core  storage  unit. 

Suppose  we  theoretically  require  a  sample  of  10,000  observations. 
For  y  •  0.5  this  would  result  in  the  collection  of  11,746  observations, 
and  for  y  ■  0.3333,  10,195  observations  would  be  collected.  Dividing 
by  the  processing  times  per  observation,  we  note  that  42.9  and  54.2 
seconds  are  consumed  when  Y  -0.5  and  Y  •  0.3333,  respectively.  Com¬ 
pared  to  the  time  cc..aumed  in  most  simulation  experiments,  this  dif¬ 
ference  for  the  two  values  of  y  is  negligible. 

The  fourth  example  is  a  single-server  queueing  problem  with  in¬ 
dependent  and  exponentially  distributed  interarrival  and  service  times 
aad  a  first-come,  first-served  queueing  discipline.  The  mean  inter- 
arriv„  and  service  times  are  Xj  -  0.?5  and  *  C.225,  respectively, 
so  that  the  activity  level  is 


(7.11) 


p  «  -  0.9  . 
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From  theory  we  know  that  the  mean  queue  length  is 

(7.12)  n  =  p/(l  -  p)  =  9  , 

and  its  variance  is 

(7.13)  RQ  -  p/(l  -  p)2  -  90  . 

Using  the  results  of  [5  and  12],  the  variance  of  the  sample  mean  de¬ 
rived  from  observing  queue  length  for  a  time  interval  (0,  T)  is  for 
large  T 

(7.14)  Var  (3^)  ~p(l  +  p^/fU  -  p)4T]  »  3847. 5/T  , 

where  T  is  measured  in  the  same  time  units  as  A,  and  A,,. 

We  recorded  queue  length  at  N  ®  T  unit  intervals.  Since  some 
variation  is  eliminated  by  this  discrete  sampling,  we  may  regard 
3847. 5/N  as  an  upper  bound  on  the  variance  of  X„ ,  the  resulting  sam- 
pie  mean. 

As  a  reliability  criterion  we  chose 

(7.15)  Pr  (jx^  -  u(  <  cu)  1  -  8  , 

where 

(7.16)  c  -  .3  , 

B  -  .10  . 

Then  using  (4.5),  the  upper  bound  on  N*  was  1420.  Three  sets  of  rep¬ 


lications  were  run,  each  with 
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(7.17)  a  -  0.025  , 

X  -  0.3333  , 

R  -  5  . 

R  was  chosen  to  emphasize  that  It  need  not  be  fixed  at  4  as  in  the 
previous  examples.  The  initial  sample  size  in  the  first  set  was 
M  ■  250,  in  the  second  M  *  500,  and  in  the  third  M  *  1000. 

The  relevant  sample  averages  are  presented  in  Table  2.  Two  ten¬ 
dencies  are  noteworthy.  The  first  is  the  constant  underestimation  of 
the  confidence  interval  width.  This  is  no  doubt  partially  due  to  the 
use  of  sampling  theory  appropriate  for  an  absolute  reliability  criter¬ 
ion  with  a  relative  criterion.  The  second  point  is  the  improved  es¬ 
timate  of  sample  size  as  the  initial  sample  size  gets  larger.  Clearly 
the  choice  of  initial  sample  size  makes  a  difference  in  the  final  re¬ 
sult. 


Table  2 

AVERAGES  FOR  QUEUEING  PROBLEM 


Theore tical 

Sample 

M=250 

M-500 

M-1000 

Mean 

9 

8.21 

:  .15 

8.83 

Lower  confidence  point 

6.3 

5.86 

5.89 

6.41 

Upper  confidence  point 

11.7 

10.56 

10.42 

11.22 

Interval  width 

5.4 

4.70 

4.53 

4.81 

Variance  of  sample  mean 

2.71 

2.10 

1.90 

2.37 

Final  sample  size 

1462a 

1068 

1299 

1473 

Equivalent  degrees  of  freedom 

32 

28 

33 

30 

Bias  adjustment 

43 

32 

37 

46 

Includes  bias  adjustment. 
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To  measure  the  overall  adequacy  of  the  suggested  procedure,  we 
studied  the  number  of  times,  that  he  generated  confidence  intervals 
included  the  true  mean  of  9.  These  results  are  presented  in  Table  3, 
where  they  are  dichotomized  by  replications  that  terminated  after  at 
roost  one  iteration  and  by  those  that  did  not.  That  is,  if  the  experi¬ 
ment  terminated  on  the  first  statistical  analysis,  then  no  iterations 
occurred . 


Table  3 


CONFIDENCE  INTERVALS  FOR  QUEUEING  PROBLEM 


f  -  - - 

Include  Mean 

Sample 

Divided  by 

M 

Iterations 

Replications 

Sample 

Expected 

Expected 

250 

mam 

29 

8 

26.1 

0.31 

H3H 

71 

58 

63.9 

0.81' 

wEsm 

100 

66 

90  9 

0.73 

500 

<  i 

36 

15 

mm a 

>  i 

64 

56 

|  0.97 

Total 

100 

71 

90.0 

1,000 

mam 

54 

36 

48.6 

BP 

46 

43 

41.4 

1.04 

i 

100 

79 

_ 1 

90.0 

250 

i  Modified 

28 

o 

72 

mm 

mSSrn 

1  Total 

100 

— 

i  90.0 

■H Hi 

Since  8  “  0.1,  we  expect  the  confidence  interval  to  include  the 
mean  in  90  percent  of  the  replications.  The  last  column  shows  the 
ratio  of  actual  to  expected  inclusions.  Notice  that  for  zero  or  one 
iteration  the  results  are  considerably  poorer  than  for  greater  than 
one  iteration.  As  the  initial  sample  size  increases,  there  is  clear 
improvement  in  the  results  for  zero  or  one  iteration  as  well  as  for 
the  remaining  category. 
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A  priori  we  generally  have  little  knowledge  of  what  an  appropri¬ 
ate  sample  size  is  for  the  experiment  to  be  undertaken.  Therefore  we 
should  be  suspicious  of  results  generated  by  zero  or  one  iteration 
since  early  termination  implies  that  our  initial  guess  is  close  to  the 
correct  answer,  a  ramcte  possibility. 

One  way  to  guard  against  premature  stopping  is  to  require  a  mini¬ 
mum  of  two  iterations  per  experiment.  To  test  the  effect  of  this  con¬ 
straint,  we  performed  the  one-hundred  replications  again  with  H  -  250, 
with  the  added  requirement  that  if  any  experiment  terminated  with 
N  -  K  useful  observations  on  the  zeroth  or  first  iteration  ws  added 
N  -  K  to  the  required  sample  size  and  continued  the  experiment.  The 
results  are  shown  in  the  last  three  rows  of  Table  3  and  axe  to  be  com¬ 
pared  with  the  corresponding  first  three  rows.  Note  the  significant 
improvement  in  experiments  that  formerly  required  less  than  two  iter¬ 
ations.  It  is  clear  that  an  ad  hoc  rule  such  as  the  one  just  described 
improves  the  statistical  properties  of  the  experimental  results. 

As  Table  3  shows,  the  percentage  of  confidence  intervals  that 
include  the  mean  generally  falls  below  that  expected  from  theory.  A- 
part  from  premature  stopping,  we  attribute  this  difference  partially 
to  the  use  of  a  relative  criterion  with  a  distribution  theory  for  ab¬ 
solute  criterion,  and  to  the  approximation  involved  in  using  theory 
developed  for  independent  observations  for  au toco r re la ted  ones. 


One  remaining  and  no  doubt  crucial  consideration  is  the  assump¬ 
tion  of  normality  of  X^.  *e  have  used  it  because  it  is  intuitively 
plausible  and  convenient.  It  is  therefore  instructive  to  study  the 
sample  size  problem  without  this  assumption  to  understand  the  advan¬ 
tages  of  the  normality  assumption.  Using  Chebyshev’s  Inequality,  we 
have  in  the  present  setting 

(7.18)  Pr  (|jL  -  uj  ±  kvWn)  -  1  -  1/k2  . 


For  the  queueing  problem  we  have 

cu  *  k»4i/n  *  2.7  , 

(7.19)  m  *  3847.5  , 

N*  =  1420  , 


so  that 


1/k2  -  - —  -  0.3704  , 

(2.7)  N* 


(7.20) 


Pr  (1^*  -  ui  <  2.7)  >  0.6296  . 


This  means  that,  regardless  of  the  distribution  of  the  probabil¬ 

ity  is  at  least  0.6296  that  deviations  from  the  mean  are  no  greater 
than  ±2.7.  Our  results  In  Table  3  are  clearly  better  than  this  low 


boundary. 
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Suppose  we  wish  to  make  no  distributional  assumption  about 
and  to  use  Chebyshev's  Inequality  to  determine  the  required  sample 
sire  N*  for  1  -  1/k2  -  0.90.  Then 


k  «  10  , 


so  that 


Pr  (jx^  -  uj  <  2.7)  i  0.9  , 

(7.21)  k2a/N*  -  (2. 7) 2  , 

K*  -  -  5278  . 

(2.7)Z  (2.7 r 

Therefore,  **e  require  an  additional  (5278-1420)  =  3858  observations 
to  meet  the  probability  requirement.  In  general,  the  use  of  Chebyshev’s 
Inequality  is  a  luxury  that  most  users  prefer  to  forego,  given  the  in¬ 
tuitive  plausibility  of  the  normality  assumption. 
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8.  THE  MINIMUM  VARIANCE  ESTIMATOR 


In  Sec.  6,  a  knowledge  of  the  autoregressive  coefficients 

{bs;  s  =0,  ....  pj  enabled  us  to  reduce  the  bias  in  the  sample  mean 

X  .  It  is  only  natural  to  inquire  whether  or  not  a  knowledge  of  these 
N 

coefficients  can  enable  us  to  derive  an  estimator  of  the  population 
mean  u  that  has  smaller  ariance  than  has  for  a  given  sample  size 
N.  The  answer  is  yes  and  it  is,  in  fact,  possible  to  derive  the  mini¬ 
mum  variance  estimator,  which  we  do  here.  To  develop  the  idea,  we 
neglect  the  bias  due  to  initial  conditions  and  correct  for  it  later. 
Consider  the  mean  estimator 


(8.1) 


**  -£*a  - 


such  that 


(8.2) 


XX 


E(V  -  *  , 


(8.3) 


v«<v  ■  £  w..t- 


s .  t=l 
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tfe  define  the  lxN  vector  of  weights. 


(8.4a) 


ij,  =  o,.  •••»  V 


and  the  NxN  autocovariance  matrix. 


(8.4b) 


I  = 


R0  *1 


R1  Ro 


V 


*}*-] 


R» 


so  that 


(8.4c) 


v«<V  - 1*  S  si  ■ 


An  array  of  the  form  (8.4b)  is  called  a  Toeplitz  matrix  and  has  a 
number  of  desirable  properties,  one  of  which  will  benefit  us  shortly. 

We  desire  to  choose  6„  to  minimize  the  variance  (8.4c)  subject 
to  the  condition  (8.2)  .  The  weights  are  then 


(8.5) 


6  - 

-N 


(h  ’ 


where  X„  i3  a  lxN  vector  of  ones.  Then 


(8.6) 


V“<S>  - l/(^  ■ 
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The  weights  are  tunetions  oi  the  inverse  of  E.,  and,  at  first 
glance,  one  night  expect  this  array  to  be  beyond  us  unless  we  know 
all  S  aucocovariances  in  (8 -4b).  Luckily,  a  knowledge  of  the  auto¬ 
regressive  coefficients  ib  ;  s  =0,  ...»  pi  and  the  residual  variance 

'  s 

2  ■!  ik 

j  suffices  to  derive  E,.  .  Let  rJ  be  the  e lessen t  of  the  inverse  in 

row  j  +  1  and  colunn  k  +  1.  Using  a  result  of  Whittle  [12,  p.73j, 

one  nay  show  that 


aik  =  (i/=2)  -T1, 


(bj-rbfc-r  •  Wjbr«-k>  J’  k  *  “•  - 


S  -  1  , 


(8.7) 


b^  =  0  r  <-  0,  r  >  p. 


Then  the  weights  are 


r 

-  5  .  .  -  ^  '  b  /[  S  ' (N-2s)b  i  k  =*  0 ,  —  ,  p  -  1 
N-k  /  j  s  s' 


s=0  s=0 


(8.8) 


- b/-  *  2s)b/j 


k— p,  • • • t  ^  “  P  > 


assuming  chat  N  >  2p  +  i.  The  variance  is  Chen 


r 

(8.9)  VarCXjP  =  a2/[b^  (N  -  2s)bsJ 


Ac  lease  three  points  deserve  mention.  First  we  note  that 


(8.10) 


^  N  var(XN,)  -  a~Jb2  , 


so  that  Xj.,  ;.he  cos  ventional  estimator,  is  asymptotically  minimum  var¬ 
iance.  More  important,  however  is  Che  second  observation  that  (8.9) 
is  an  exact  formula  for  the  variance  in  contrast  to  using  the  limit 
~2/( Nb2)  for  Xj.-  The  third  point  is  that  for  a  scheme  of  order  p,  observa¬ 
tions  p,  ...,  S  -  p  receive  eaual  weight  so  chat  the  effect  of  autocor¬ 
relation  is  felt  at  the  ends  of  ""he  time  series. 

To  compare  Xj.  and  a  furthe”  example  is  helpful-  Consider  the 
first-order  autoregressive  scheme  (2.30)  again.  We  have 


3l  --  SN  =  l/[M(i  -  g)  2gj  , 

(8.11)  SR  -  (1  -  g)/[N(l  -  g)  +  2gj  k  -  2,  ....  H  -  1  , 

Var^J  -  1/[H(1  -  g)2  +  2g(l  -  g)j  . 

Notice  that  for  0  <  g  <  1  the  first  and  last  observations  receiv  •  the 
largest  weights.  Since  the  first  observation  is  subject  tc  the 
greatest  bias  frea  Che  initial  conditions,  it  is  important  to  remove 
this  bias  when  the  minimua  variance  estimator  is  used.  By  contrast, 
we  note  that  for  -1  <  g  <  0,  the  first  and  last  observations  receive 
less  weight  than  the  remaining  observations. 

Let  N  be  the  required  sample  sire  for  a  specified  reliability 
using  (2.5a),  and  let  be  the  required  sample  size  using  (8.1) 


Fig.  4 — Comparison  of  and  3^  for  first-order  autoregressive  schemes 


Figure  4  shows  N/N  for  values  of  N,  and  g  *  0.95,  0.99.  For  g  *  0.95, 
<v. 

the  advantage  of  X^  is  most  outstanding  for  sample  sizes  of  less  than 

v 

100.  The  desirability  of  therefore  is  most  apparent  when  high  unit 
collection  costs  require  a  limit  on  the  sample  size.  Perhaps  more  im¬ 
portant  is  the  fact  that  the  expression  for  the  variance  is  known  ex¬ 
actly  for  any  order  scheme. 

The  bias  adjustment  can  be  introduced  here  simply  by  discarding 
the  first  K  observations  and  using  (8.1),  with  N  replaced  everywhere 
by  N  -  K,  provided  that  N  -  K  i  2p  +  1.  The  new  mean  estimator  is 
then 


(8.12) 


Y> 


N-K 

■2>a* 

t=i 


» 


with  variance 


P 

Var(XN  R)  -  aZ/[b^(N  -  K  -  2s) bg]  . 


(8.13) 
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"§C 

Suppose  that  it  is  desired  to  compute  N  satisfying  the  absolute 
reliability  criterion  (4.1) -  Using  (7.9)  we  have 

P 

(8.14)  N*  -  m(Q/c)2  +  (2/b)  ^sbg  . 

s-0 

If  the  relative  criterion  (4.9)  is  used,  we  have 

P 

(8.15)  H*  -  o[Q/(cp)32  +  (2/b)7>a  . 

s*0 

If  the  only  specification  is  that  the  variances  be  less  than  or  equal 
to  V ,  then 

p~ 

(8.16)  N*  -  n/V  +  (2/b)  ^sbg  . 

s-C 

This  is  the  criterion  shown  in  the  flow  chart  in  Fig.  5  with  the  bias 
correction.  It  will  be  helpful  to  note,  when  computing  the  6^'s,  that 

SK+1  "  9N  “  e’  * 

(8.17)  6S+l+k  "  eN-ic  "  ®*fk  +  9’bk  k-1,  ...,p-l, 

\+l+k  “  ®K+p  k  "  p»  •**•  N  *  P  * 

P 

0'  -  1/C  ^P(N  -  2s) bs]  . 

«■ 0 


I 


f: 
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9.  CONCLUSIONS 

This  Memorandum  has  presented  a  number  of  ideaa  which,  when  taken 
together,  enable  the  user  to  derive  results  from  his  simulation  experi¬ 
ment  with  a  specified  accuracy.  Moreover,  the  results  are  obtainable 
without  user  interactions  with  the  ongoing  experiment. 

Zn  contrast  to  the  spectral  approach  [5,6],  which  often  requires 
the  computation  of  a  large  number  of  autocovariances  to  estimate  the 
variance  of  the  sample  mean,  the  suggested  autoregressive  approach  re¬ 
quires  R  sample  autocovariances  where  R  need  not  exceed  4  or  5.  The 
autoregressive  approach  also  supplies  more  objective  estimates  of  the 
equivalent  degrees  of  freedom  that  in  turn  pemlt  the  computation  of 
confidence  Intervals  explicitly  acknowledging  ths.  estimated  variance 
of  the  same  mean.  The  ability  to  correct  for  bias  due  to  initial  con¬ 
ditions  is  also  noteworthy. 

The  Memorandum  has  dealt  with  the  problem  of  determining  size  for 
estimating  the  mean.  An  ancillary  problem  that  deserves  attention  is 
how  to  estimate  the  mean  efficiently.  To  this  end.  Sec.  8  has  described 
the  minimum  variance  unbiased  estimator,  which  turns  out  to  be  a  frac¬ 
tion  of  the  autoregressive  coefficients.  No  doubt  such  results  will 
interest  tho» *-  concerned  with  variance  reduction  techniques. 

The  flow  charts  provide  procedural  assistance  for  implementing 
the  technique  described.  The  omission  of  an  example  using  the  minimum 
variance  estimator  is  deliberate,  since  we  prefer  to  learn  more  about 
the  resulting  sampling  properties  before  using  it  extensively. 


-54- 


REFERENCES 


I.  Cox,  D.  R. ,  and  P.A.W.  Levis,  The  Statistical  Analysis  of  Series 

of  Events,  Methuen,  London,  1966. 

x.  Cox,  D.  R. ,  and  H.  D.  Miller,  The  Theory  of  Stochastic  Processea, 
John  Wiley  and  Sons,  New  York,  1965. 

3.  Durbin,  J.,  "The  Fitting  of  Time-Series  Models,"  Rev.  Int.  Inst. 

Slat.,  Vol.  28,  1960,  pp.  233-244. 

4 .  Fishman ,  G .  S . ,  Problems  in  the  Statistical  Analysis  of  Simulation 

Experiments:  The  Comparison  of  Means  and  the  Length  of  Sample 
Records ,  The  RAND  Corporation,  RM-4880-PR,  February  1966.  Also 
publish'!  in  Comm.  ACM,  Vol.  10,  February  1967,  pp.  94-99. 

5.  _ ,  and  P.  J.  Kiviat,  "The  Analysis  of  Simulation  Generated 

Time  Series."  Mgt.  Sci.,  Vol.  13,  No.  7,  March  1967,  pp.  525-557. 

6.  Fishman,  G.  S.,  "The  Allocation  of  Computer  Time  in  Comparing  Simu¬ 

lation  Experiments,"  Opns.  Res.,  Vol.  16,  No.  2,  March-April  1968 
PP.  280-295. 

7.  Fraser,  D. ,  Nonparrmetric  Methods  in  Statistics,  John  Wiley  and 

Sons,  New  York,  if 57. 

8.  Jenkins,  G.  M. ,  ani  D.  G.  Watts,  Spectral  Analysis  and  Its  Appli¬ 

cations  ,  Holden  Daj ,  San  Francisco,  1968. 

9.  Kendall,  M.  G.,  and  A.  Stuart,  The  Advanced  Theory  of  Statistics. 

Vol.  3,  Hafner,  London,  1966. 

10.  Kiviat,  P.  J.,  R.  Villanueva  and  H.  M.  Markowitz,  The  SIMSCRIPT  II 
Programing  Language,  The  RAND  Corporation,  R-460-PR,  October 
1968. 

II.  Malinvaud,  E.,  Statist ica.  Methods  in  Econometrics,  Rand  McNally, 

Chicago,  1966. 

12.  Morse,  P.,  "Stochastic  Properties  of  Waiting  Lines,"  Opns.  Res., 

Vol.  3,  No.  3,  August  1955,  pp.  255-261. 

13.  Quenouille,  M.  H.,  "A  Large-Sample  Test  for  the  Goodness  of  Fit  of 

Autoregressive  Schemes,"  J.  Roy.  Stat.  Soc.,  Series  A,  Vol.  110, 
1947,  pp.  123-219. 

14.  U.S.  Department  of  Commerce,  Handbook,  of  Mathematical  Functions, 

National  Bureau  of  Standards,  Washington,  1964. 

15.  Whittle,  P.,  Prediction  and  Regulation  by  Linear  Least -Square  Meth¬ 
ods,  The  English  Universities  Press,  London,  1963. 

Wold,  H.,  A  Study  in  the  Analysis  of  Stationary  Time  Series,  2d  ed. 
Almquist  and  Wicksell,  Stockholm,  1954. 


16. 


DOCUMENT  CONTROL  DAT.'' 


[I  ORIGINATING  ACTIVITY 


THE  RAND  CORPORATION 


2a  REPORT  SECURITY  CLASSIFICATION 

_ UNCLASSIFIED _ 

26.  GROUP 


RE»OPT  T.TLE 

DIGITAL  COMPUTER  SIMULATION’:  ESTIMATING  SAMPLE  SIZE 

AUTHOR(S)  {Loti  namt,  first  noma, initial) 

Fishman,  George  S- 


REPORT  DATE 

August  1969 
CONTRACT  OR  GRANT  No. 

F44620-67-C-0045 

c  AVAILABILITY/ LIMITATION  NOTICES 


6o.  TOTAL  Ho.  OF  PAGES 

_ 61 _ 

8.  ORIGINATOR  S  REPORT  No 


RM-5S66-PR _ 

1 9b.  SPONSORING  AGENCY 


iSb.No.  OF  REFS. 


DDC-1 


J.  ASSTRACT 


An  algorithn  for  automatically  estimating 
and  collecting  the  sample  size  required 
for  statistical  precision  in  a  computer 
simulation  experiment  while  the  simula¬ 
tion  is  running.  The  algorithm,  which 
would  be  incorporated  directly  into  the 
computer  routines,  w’ould  relieve  an  in¬ 
vestigator  of  the  burden  of  first  esti¬ 
mating  the  variance  of  the  sample  mean 
from  a  data  sample  obtained  from  a  trial 
run,  then  estimating  the  sample  size  ne¬ 
cessary  for  the  specified  confidence  in¬ 
terval,  and  finally  collecting  that  many 
more  observations  in  a  successive  simula¬ 
tion  run.  The  underlying  probability 
model  is  autoregressive:  it  -would  depend 
on  an  autoregressive  representation  of  the 
sequence  that  considers  each  observation 
as  a  linear  combination  of  past  observa¬ 
tions  plus  an  unccr related  random  resid¬ 
ual.  This  approach  need  not  require  more 
than  4  or  5  autocovariance  computations 
to  estimate  the  variance  of  the  sample 
mean.  A  flowchart  is  included  to  aid  in 
building  the  technique  into  simulation 
programs . 


United  States  Air  Force 
Project  RAND 

II  KEY  WORDS 


Computer  simulation 
Statistical  methods  3nd  processes 


